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A variational formulation for tiie calculation of interacting fermion systems based on the density- 
matrix functional theory is presented. Our formalism provides for a natural integration of explicit 
many-particle effects into standard density-functional-theory based calculations and it avoids am- 
biguities of double-counting terms inherent to other approaches. Like the dynamical mean-field 
theory, we employ a local approximation for explicit correlations. Aiming at the ground state only, 
trade some of the complexity of Green's function based many-particle methods against efficiency. 
Using short Hubbard chains as test systems we demonstrate that the method captures ground state 
properties, such as left-right-correlation, beyond those accessible by mean-field theories. 
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I. INTRODUCTION 

One of the most successful theories in solid-state 
physics is the density-functional theory (DFT)ii^. Its 
present implementations predict properties of materials 
with an accuracy, which makes it a standard tool for 
material scientists both in fundamental research and in 
industry. 

One of the major deficiencies of currently available 
density functionals is the description of materials with 
strong electron correlations. Strong electron correlations 
develop when the interaction between the electrons dom- 
inates over the kinetic energy. While the kinetic energy 
favors delocalized electron states and well-defined energy 
bands in reciprocal space, strong interactions may lead 
to a complete breakdown of the band picture. 

Strongly correlated materials often contain 3d 
transition-metal elements or rare-earth elements with 
partially filled f-electron shells. The d- and f-electrons 
are on the one hand localized core-like states, but, on 
the other hand, they are located energetically in the re- 
gion of valence states. Thus, their nature is right at the 
border between localized and delocalized behavior. 

Strong correlations are responsible for a wealth of 
new properties with potential technological relevance. 
Among them are hcavy-fermion systems^, Mott insula- 
tors and the Mott-Hubbard metal-insulator transition 
(MHMIT)ii. high-temperature superconductors^ and the 
colossal magnetoresistance^. More recently, the vicinity 
to phase transitions in those systems has opened the field 
of quantum criticality^. 

It is a fascinating challenge to unravel the underlying 
mechanisms behind these phenomena and eventually to 
rationally design applications exploiting them. A theo- 



retical treatment of these systems, however, turns out to 
be a rather big challenge and has been at the frontier of 
research in condensed matter physics since nearly half a 
century.— The use of explicit many-body tools is nearly 
impossible for complex materials like perovskites. There- 
fore, the many-body community tries to identify the min- 
imal set of degrees of freedom relevant for the physical 
properties and sets up a model hamiltonian for these de- 
grees of freedom. The best known of these models is the 
Hubbard model^"— , which has been the working horse 
for correlation effects in transition metal compounds for 
almost 50 years. 

A first attempt to combine DFT with techniques from 
many-body theory was made with the GW methodic, 
which supplements DFT with certain classes of Feynman 
diagrams. While the GW scheme was applied to semi- 
conductors with some successii, it still fails to describe 
the physics of the MHMIT. 

A partial success along this line was achieved by the 
so-called LDA-I-U approach, which complements the den- 
sity functionals by Hartree-Fock type local correlationsi^. 
Hybrid density functionals^^, which replace part of the 
exchange energy with the exact non-local exchange en- 
ergy, capture essentially the same physics. These func- 
tionals lead to a large improvement of the band gaps, 
and transition-metal oxides appear correctly as antiferro- 
magnetic charge-transfer insulators, where conventional 
GGA-type functionals erroneously predict metals or a 
small d-d band gapjii^ It is, however, important to em- 
phasize that within LDA-I-U one can only describe mate- 
rials with ordered ground states satisfactorily, while the 
MHMIT in the paramagnetic state, as in V2O3; cannot 
be explained with this method. 

A quantum leap in the understanding of the MHMIT 



was achieved by the invention of the dynamical mean- 
field theory (DMFT)^^, which provides an efficient and 
reliable tool to approximately solve the correlation prob- 
lem for Hubbard type models with local interactions^. 
DMFT keeps the local dynamics intact, which turns out 
to be essential for the MHMIT. Hence, this theory can 
actually describe the competition between the Fermi liq- 
uid and the interaction driven insulato r^^i^^ , and also 
allows for the inclusion of ordered phases beyond stan- 
dard Hartree theoryi^"— . The price to pay is the ne- 
glect of non-local physics and that DMFT fails to cap- 
ture phenomena like unconventional superconductivity or 
quantum criticality based on spatial fluctuations^. Note, 
however, that these effects can be re-included to a certain 
extent by extending the original theory^S 

Rather soon it was realized, that the DMFT can be 
combined with DFT in a fashion similar to GW<^"— In- 
deed, early successes of the method comprise the proper 
description of spectral properties of (Sr,Ca)V03, the 
MHMIT in V2O3, the energetics of the volume collapse 
in cerium^^, and the lattice properties of plutoniumM, to 
name a few. 

However, as the DMFT is model-hamiltonian based, 
it requires the definition of a Hubbard model prior to 
its application. This means, that one needs to extract 
from DFT (i) the hopping parameters of the Hubbard 
model and (ii) the local Coulomb parameters. Both are 
highly ill-defined quantities, because the results strongly 
depend on the basis set and the sub-manifold of orbitals 
taken into account. Here, a large effort has been put 
forward during the past years to define a fairly universal 
interface between the wave-function based DFT and the 
Green function based DMFT^"^. 

Besides these technical problems, there are further 
caveats of the method. Firstly, the DFT does already 
include certain correlation effects. These have to be 
"subtracted" from the DFT prior to the calculation of 
tight-binding and interaction parameters to avoid double- 
counting in the DMFT calculation. This subtle aspect is 
one of the most challenging problems in the combination 
of DFT and DMFT so far, because nobody knows to 
what extent and in what form in the diagrammatic lan- 
guage correlations effects are included in DFT and how 
to actually properly remove thero^^. 

Secondly, within a many-body formulation, the inter- 
actions within the subspacc of the correlated orbitals will 
be screened by the other electronic degrees of freedom in 
the spirit of GW— . This means, that the actual interac- 
tions become retarded. This poses a serious challenge for 
the computational tools necessary to solve the effective 
quantum impurity model arising within DMFT. 

And last, but not least, the combination of DFT 
and DMFT cannot be put rigorously on a variational 
basis. One can, in principle, formulate an extremal 
principl o^^'^^ involving both wave function for DFT and 
Green functions for DMFT, but the actual connection 
between these two worlds and a proper subtraction of 
double counting terms is not obvious. 



Therefore, a method, which is fully variational and al- 
lows to include both correlations as they are described 
by conventional density functionals and a reasonable and 
controlled approximation for strong local correlations, is 
clearly desirable. The aim of this paper is to propose 
such an approach. The guiding idea is to reformulate the 
variational problem in terms of the one-particle density 
matrix, which is known as rcduced-dcnsity-matrix func- 
tional theory (rDMFT). This problem is of the same com- 
plexity as a complete many-particle calculation. There- 
fore we introduce a local approximation to avoid the di- 
mensional bottleneck. 

The origin of rDMFT goes back to Gilbert's theorent^. 
Gilbert proved that every ground-state observable can be 
written as a functional of the one-particle density ma- 
trix. Analogously to the development in the field of 
density-functional theory, research is directed towards 
the formulation of approximate density-matrix func- 
tionals. Most applications use a variation of Miillcr's 
functional^^Ti^. Recent functionals have been able to 
describe the Mott band gap of transition metal oxides 
within a non-magnetic description^!, while density func- 
tionals tend to form antiferromagnets instead. This prob- 
lem is related to the broken-symmetry states of DFT. 

In this paper, we explore the main principles of the 
method. We do not yet approach the question on how to 
exploit the local nature of the interaction in the Anderson 
impurity models. 

The paper is organized as follows. In the following 
section |lll we present the formulation of the theory in 
detail. In section IIIH we sketch the evaluation of the 
reduced-density-matrix functional using a dynamical op- 
timization similar to the Car-Parrinello method^. The 
workings of our theory will be explored in section IIV| 
where its results for short Hubbard chains are compared 
with rigorous many-particle calculations. 



II. THEORETICAL FOUNDATIONS 

A. Notation 

While our goal is to incorporate many-particle effects 
into a general electronic structure method such as the 
projector augmented-wave (PAW) method method^^, the 
first step will be to map the electronic structure onto a 
basisset of local spin orbitals Xai^.a) = (f,a\xa)- The 
spin coordinate a £ {f, i} can assume the values spin-up 
and spin-down. Each orbital is a two-component spinor, 
even though we may choose the orbitals such that they 
only have one non-zero spin component. 

In order to decompose an arbitrary one-particle wave 
function j-f/;) into the contribution of the orbitals \xa) 
according to 



W ='^\Xa){Tra\llj) 



(1) 



we introduce so-called projector states ItTq). The projec- 
tor states TTair^a) = (r, (tItTq) obey the bi-orthogonality 
condition 



has the matrix elements 



(tTqIX/j) =Sa,l3- 



(2) 



The identity Eq. [T] is valid whenever the wave function 
can be expanded completely into the orbital set. For a 
complete orthonormal basisset of local orbitals, the pro- 
jector states can also be chosen identical to the local or- 
bitals. However, in the most general form, the projector 
functions differ from the orbitals. 

The projector states ItTq.) are conceptually related to 
the projector functions of the PAW method^. They dif- 
fer from those used for the augmentation, because they 
provide the amplitudes of the local orbitals instead of 
partial waves amplitudes provided by the latter. How- 
ever, in the PAW method the projector states ItTq.) for 
the local orbitals are naturally constructed as superposi- 
tion of the projector states for the partial waves. 

The creation and annihilation operators in terms of lo- 
cal orbitals can be expressed by the field operators ijj{r, (j) 
and ^^ (r , a) via 



4 = X! / d^ril)\f,a)TTair,a) 

a •' 
Ca = ^ / (i^r 7r*(f, cr)Vi(f, cr) . 



The back transform is 

a 



(3) 
(4) 

(5) 
(6) 



Creation and annihilation operators obey the anti- 
commutator relation 



(TTalTT/j) 


(7) 





(8) 


0. 


(9) 



[CLC/3]. 



-t 't 



[Ca,C/3]_ 



Note, that the anti-commutator between creators and 
annihilators, i.e. the overlap between projector states, is 
not necessarily equal to Sij. This is due to our general- 
ization to non-orthonormal local orbitals. 

The Hamilton operator is given by 



H = h + W . 



(10) 



where h is the one-particle part of the Hamilton operator, 
containing kinetic energy and external potential. W is 
the Coulomb repulsion between the electrons. 
The non-interacting Hamiltonian 



a, 13 



(11) 



ha,p = {Xa\-^ + VextlXli) ■ 



(12) 



The interaction energy has the form 

^ = 2 ^ W^a,/3,7,5CJ,C]5C5C^ (13) 

with the matrix elements 

e^xUr,'^)x*fiir',(T')X'rir,cr)xs{r',(T') 
X .(14) 



Note the order of the indices of Wa^p.^^s in Eqs. [T3] and 

m 



B. Reduced-density-matrix- functional theory 

As shown by Gilbert^, every ground-state property of 
an electron gas can be expressed as functional of a one- 
particle reduced density matrix. Let us review the basis 
of the theory: 

The ground-state energy for a given external potential 
can be written as 

£;(h) = min ( ($1 J2 K.,pcicp + W\<^) 



£{{m)-l\ 



(15) 



The minimization is performed over all fermionic many- 
particle wave functions |<I>). £ is the Lagrange multiplier 
for the normalization constraint of the wave function. 
The particle number is not restricted and the Fermi level 
is placed at the energy zero: that is, we implicitly employ 
a grand-canonical description for the electrons: Electrons 
are taken from some vacuum level that defines our energy 
zero and added to the system as long as this addition 
lowers the total energy. 

Let us perform a Legendre transform of the total en- 
ergy, Eg. llSi with respect to the one-particle Hamiltonian 



F^{p) = min [£:(h) - Tr(hp)] 



(16) 



The minimum condition of Eq. [16] with respect to h, 
dE 



Pa,/S(h) 



dhp.a 



$|clca|$) 



^/3 



(17) 



identifies the conjugate variable p with the one-particle 
reduced density matrix of the ground state |<f>) of the 
Hamiltonian h + W. 



The result of the Legendre transform is the rcduccd- 



density-matrix functional 



F 



w 



(p) = min 

|*),h,£ 



($|iy|$) 



aJ3 



,/3((*l4c/3|$)-p^„ 



£(($|$)-1) 



(18) 



The reduced-density- matrix functional _F [p] is univer- 
sal in the sense that it is an intrinsic property of the inter- 
acting electron gas and that it is independent of the ex- 
ternal potential. We introduced the interaction W as su- 
perscript in the notation for the reduced-density-matrix 
functional because we will use different interactions in 
the course of this paper. 

The one-particle Hamiltonian h required in Eq. [16] 
to produce a specified density matrix p is obtained, up 
to the different sign, as derivative of the density-matrix 
functional with respect to the density matrix, namely as 



/Ja,/3(P) = - 



dF 



w 



dpp.a 



(19) 



The range of allowed arguments p for the universal 
functional is only a subset of all hermitian matrices in 
the one-particle Hilbert space. The range is limited to 
the matrices that can be constructed by Eq. [17] with any 
fermionic many-particle state 1$) from the Fock space. 
Matrices of this form are called N-representable^ 

According to a theorem due to Coleman^, the eigen- 
values of all N-representable matrices lie between zero 
and one and all matrices with eigenvalues between zero 
and one are N-representable. 

The eigenstates ji/'n) of the one-particle reduced den- 
sity matrix are called natural orbitals^ and the eigen- 
values /„ are their occupations. The one-particle density 
matrix expressed in terms of natural orbitals and occu- 
pations is given by 



Pa, 13 = y^(7ra|^n)/n(-0n|7rff) 



(20) 



Hence, the total energy can be expressed by a Legendre 
back-transform as 



£:(h) = min 

{|V-„),:.„,A„,„} 



- ^ A„_„i ((V'rJV'm) - 5,n,n) ■ (21) 



In order to respect the limited range of values for the 
occupations we introduce the transformation f{x) — 
sin (^x) from an unrestricted real variable x to the in- 
terval [0, 1]. The An,m are the Lagrange multipliers for 
the orthonormality constraint of the natural orbitals. 



Starting from the basic formulation of reduced-dcnsity- 
matrix-functional theory, the usual route taken is to find 
and explore approximate but analytic expressions for the 
density-matrix functional. We follow a different route, 
namely to determine the functional on the fly using the 
explicit constrained-search method of Levy^. 



C. Local approximation 

Up to now, the formulation has been exact. How- 
ever, the evaluation of the density-matrix functional us- 
ing the constrained search algorithmic requires the solu- 
tion of a full many-particle problem. Because the many- 
particle Hilbert space grows exponentially with system 
size, we introduce approximations that avoid this dimen- 
sional bottleneck. 

To this end, we introduce clusters Cr of orbitals for 
which the correlations will be treated explicitly. Such a 
cluster may include all orbitals from the d- or f-clcctron 
shell of an atom; alternatively, it may include all orbitals 
centered on a single atom, or it may include all or a 
certain subset of orbitals from several atoms. In analogy 
to the Anderson impurity model, we will refer to these 
clusters in the following as " impurities" . 

Besides the full interaction of Eq. [TS] we introduce a 
local interaction Wioc- Two electrons interact with the 
local interaction only if both reside on the same impurity, 
i.e. in the same cluster Cr of correlated orbitals. 



W-loc = 51 ^« 



Wr 



9 Yl ^'^^Pn.s clclcsc^ 



(22) 



(23) 



a,fi,f,SeCR 



Such an interaction is well-known in model-based studies 
of interaction effects. For example, the Hubbard model 
uses an interaction of this form4"— 

In order to establish the link between DFT and DMFT, 
let us also introduce an approximate density-matrix- 
functional F^prp[p] calculated from a density functional 
as 



F 



DFT 



[P] 



d-^r / d'^r' 



, e'^n(r)n{r') 
47reo|^— ''' 



E^,[n{r)] ,(24) 



where n{r) := Y.a,i3 Pa.0X}i['r)xa[r) is the electron den- 
sity. 

With the approximate functional from Eq. [M] the 
density-matrix functional with the full interaction can 



be written in the form F^ = F, 



DFT 



{F 



w 



^DFTJJ 



i.e. as DFT plus a correction describing many-particle 
effects explicitly. Irrespective of the choice of the density 
functional in Eq. [24] this form of the functional is exact, 
because it simply adds and subtracts the same term to 
and from the full functional F^ . 

The first approximation of our theory is to replace the 
interaction in the correction term with the sum of local 



interactions 



F 



w 



FL'ft+(f^-'^''-fE%^^). (25) 



pW 



On this level, the theory describes interactions on the 
impurities with the full functional, while the non-local 
parts of the Coulomb interaction are captured on the 
DFT level. 

The second approximation of our theory is to replace 
the correction term in Eq. [25]by a sum of local correction 
terms 



t ~ i'DFT 



y^ fpVVn. 



pWn 
^ DFT 



(26) 



In this approximation, the environment of each cluster is 



J 



E[V,^uN] 






p 



replaced by a non-interacting electron gas. Each term 
F'^'^ has precisely the structure of a so-called multi- 
orbital single-impurity Anderson model4^ 

For the sake of completeness, let us quote the error 
term for the approximation of Eqs. 1251 and I26 ( 



AF'' 



pW 






F, 



w 

DFT 



/ ,^DFT 



(27) 



Ve, 



2ra^ 

{'^a\^n)I{Xn){'>Pn\T^I3 



With Eq. [26] we obtain, after including the particle- 
number constraint, the following expression for the total 
energy. 






47reo|r — r'\ 
2 J 47reo|r — r'\ 



EZ'^ln^ 



(.r)]} 



n , m n ) 



(28) 



where, with /„ = /(a;„). 



D. Exchange-correlation w^ith local interaction 



and 



^(^) = I]I](^,^l^n)/„(V'„|r,ff), (29) 

(T n 

X^(7ra|-0„)/„(-0n|7r/3)(x/?|^CT) , (30) 



"flW = E E (^''^l^") 
o- a,P£CR 

>^^{Tra\^Pn)f7i{-ipn\Trp){Xls\r,Cr) ■ (31) 



Finding an expression of the exchange-correlation func- 
tional for the interaction is not straightforward. 

The exchange-correlation energy can be expressed 
by the interaction-strength averaged hole function^ 
h\{r,r') as 



Exc[n]^ d?rn{r) 



2 J 47reo|f-r'| 



.(32) 



The hole function is defined via the two-particle density 

(2) ^ 

n\ (r, r') of the ground state for a given density n{r) or 
spin density n^,a-' (f) as 



hx{r,r') 



n{r) 



n\^ (r, r') — n{r)n{r') 



(33) 



For the sake of simplicity, we write the expression for 
density functionals and not for spin-density functionals. 
The generalization is straightforward. 

The second line in Eq. [28]can be considered as a corre- 
lation correction to the density functional: In the absence 
of the correlation correction, the theory recovers the con- 
ventional expression of density-functional theory. 



The parameter A scales the interaction W. The corre- 
sponding hole function is obtained from the wave func- 
tion for this scaled interaction. The A-average accounts 
for the difference of the kinetic energy between the inter- 
acting and the non-interacting electron gas. 

Because the restricted interaction Wr is nonlocal in 
the coordinates of each particle, it needs to be modified 
before it can be used in a density functional. We propose 



to use the following model for the loeal interaetion 



.,Wr 



{r,r') 



nnir) 



n^(r') 



n^{r) 4:7reo\r-r'\ n^{r') ' 



(34) 



where n^{r) and n^{r) are defined by Egs. [501 and [5T] 
This choice for the interaction is consistent with the ex- 
pression for the Hartree energy. It can systematically 
be extended to the full exchange correlation energy by 
including the remaining pair terms. 

If we use the model for the restricted interaction from 
Eq. [31 in Eq. [31 we obtain 



^Wr 



(fr n^{f) 



- I £-' 



e'^h{r,r') n^ir') 
47reo|r — r'\ n^{r') 



(.35) 



Here, the A-averaged hole function is used. Hole func- 
tions for some of the commonly used density functionals 
have been developed^i^. 

Eq. lBSI rcflects that only electrons on the same impurity 
interact with each other. Our approximation assumes 
that the exchange-correlation hole of the Anderson im- 
purity model is identical to that of the fully interacting 
system. This is strictly valid for the exchange part, while 
it is an approximation for the correlation contribution. 

If the hole function is not available, the most simple 
approximation is 



Sfc"- d'rnl{r)e,,[nl^,{r}] 






(36) 



where Cxc is the exchange correlation energy per electron. 



E. Coulomb matrix elements 




<— > 



f oV© 



FIG. 1: Scheme to demonstrate the case where an electron 
leaves a correlated region. In the local approximation the in- 
teraction with the electron remaining on the correlated cluster 
is completely removed, explaining the need for renormalizing 
the Coulomb interaction. (Color online) 

One of the important ingredients to our theory is the 
Coulomb tensor (|23p projected onto the impurity or- 
bitals. The first idea is to calculate this tensor using the 
local orbitals \xa) with the bare Coulomb interaction. 

However, as in other local approximations, such an ap- 
proach would result in too large Coulomb parameters for 
the reason schematically depicted in Fig. [T] If an elec- 
tron leaves a correlated orbital, it is likely to be located 
nearby, so that only a portion of the interaction is lost. 



In the local approximation the Coulomb interaction is 
completely removed for a particle leaving the correlated 
impurity. Thus the energy required to add or remove an 
electron from the impurity is overestimated in the local 
approximation, respectively charge fluctuations are arti- 
ficially suppressed. Therefore, our method will require 
an appropriate renormalization of the local U-tensor to 
account for those screening channels. 

There are several possible approaches to obtain this 
effect. One can in principle enlarge the impurity and in- 
clude more extended orbitals in the correlated cluster. 
This will lead at least partially to the desired renormal- 
ization: The screening is taken over by the extended 
orbitals. Another route, which has proven to be more 
efficient, is to perform additional calculations using con- 
strained DFT— , which will also provide an estimate of 
the renormalization due to screening from other orbitals. 



F. Extension to ensemble density-matrix functional 
theory 

The formulation of the density-matrix-functional the- 
ory used so far is based on pure states, that, however, 
are not necessarily eigenstates of the particle-number 
operator. Usually, density-matrix-functional theory is 
formulated for ensembles of many-particle states. The 
generalization to ensemble rDMFT is straightforward. 
A density-matrix-functional for an ensemble of many- 
particle states has the form 

F^(p) = „min £P,(*.ri$.) 



^/i„,J^P,(<i>,|clc^|<f>,) 



PP, 



a. 13 



^ A,,, (($,!$,) -<5.,,) + a(^P,-1) 



«j 



(37) 



where Pi are probabilities for the many-particle states 
$i) . In order to ensure that the probabilities fall between 
zero and one, we represent them as square Pi = Xf of a 
real number Xi. 

Formally, the theory can be extended to finite tem- 
peratures by adding an additional entropy term —ST = 
ksT J2i Pi lii(Pi), that represents a heat bath, to the den- 
sity matrix functional. 



III. EVALUATION OF THE DENSITY-MATRIX 
FUNCTIONAL 

The reduced-density-matrix functional F^« in Eq. [28] 
is determined directly using the constrained-search for- 
malism of Levy^. The details of the methodology will 
be described elsewhere. 



In order to evaluate the density-matrix functional, 
we translate the fictitious Lagrangian methodology of 
Car and Parrinello for ab-initio molecular dynamic 



.38 



to 



many-particle wave functions |<I>). 

The problem of ab-initio molecular dynamics and of 
a many-particle problem are related, as far as they are 
applied to the search for ground states: in both cases 
one has to determine one or a few of the lowest states 
of a highly-dimensional Hamiltonian. Furthermore, the 
fictitious Lagrangian formalism lends itself to the search 
of a minimum under constraints. In the Car-Parrincllo 
method^^, the constraints are the orthonormality con- 
straints for the Kohn-Sham wave functions. In our case 
it is the norm and the requirement that the wave function 
produces a specified reduced one-particle density matrix. 

We start from a fictitious Lagrangian defined as 



^m,m 






a, 13 



1 



(38) 



In practice, each many-particle wave function is repre- 
sented by a coefficient array of the Slater determinants. 
Including a friction term, the equation of motion has 
the form 



TO$|(f>) 



-(> >a./3cJ,C/3 






W -E)\^) - |$)k.(39) 



This equation of motion is discretized using the Verlet 
algorithm^, and the constraints are taken into account 
using the method of Ryckaert et al4^. 

The Car-Parrincllo method is not only used to optimize 
the many-particle wave function for each Anderson im- 
purity model. It is also used more conventionally to opti- 
mize dynamically the natural orbitals and occupations in 



J 



the "outer loop" . This is analogous to the conventional 
Car-Parrinello method, which optimizes the Kohn-Sham 
wave functions. Both, natural orbitals and Kohn-Sham 
wave functions are single-particle wave functions. 



IV. MODEL CALCULATIONS 

In order to explore the workings of the formalism de- 
scribed in the previous sections, we performed calcula- 
tions on model systems that allow comparison to exact 
results. 

The model systems arc Hubbard chains with Ng ~ 2, 
Ns = 4:, and Ng ^ 6 sites. Their Hamiltonian is 



N, 



H = T+Y^Wr, 



(40) 



fl=i 



where T describes the non-interacting part of the Hamil- 
tonian 

T=-iYl E (4,cTCfl+i,<T + 4-Hi,<TCfl,<T) , (41) 

and Wr is the interaction on a the R-th site 

Wr = C/%,t^i?4 • (42) 

t is the hopping parameter and U is the Coulomb param- 
eter. fiR^a = Cr ^CR^a is thc particlc-numbcr operator for 
the specified orbital. 

We adapted Eq. [^S] to the Hubbard chain by ignoring 
Hartree and exchange energy terms. They cancel exactly 
on this level of the theory. 



E{N) = min ^ ^ /(x„)(V'n|r|^„) + E ^"^"(i/^"^^}) 



(43) 



where the f{xn) are the occupations defined in Eq. (|2ip . and the one-particle states \ipn) arc the natural orbitals. 
Thc density-matrix functional for for each site has thc form 






mWR^ + Y, ha.P ( ($l4c/3|*) - Pp.a )-E[ ($1$) - 1 
a, 15 



(44) 



r 



The total energy for half filling is shown in figure [5] as and how thc interaction switches off the covalent interac- 
function of thc interaction strength U. Of interest is if tion due to left-right correlation^. Left-right correlation 
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FIG. 2: Energy as function of the interaction strength of the 
2-, 4-, and 6-site Hubbard chains at half filhng. Resuhs from 
our theory (symbols) are compared with the exact result (solid 
line) and the mean- field approximation (dashed line). 



describes that two electrons forming a bond avoid being 
on the same site in order to lower their Coulomb repul- 
sion and plays a major role in chemistry, because it is vi- 
tal for understanding the dissociation of chemical bonds. 
Due to this effect, the covalent interaction energy is com- 
pletely lost as the interaction grows. To make contact to 
the correlation effects in solids, we note that left-right 
correlation describes the same physical mechanism that 
is responsible for the Mott band gap. 

For the dimer at half filling, our theory reproduces the 
exact result, because the many-particle wave function for 
the two-particle state is completely determined by the 
density matrix. For the other systems, our theory un- 
derestimates the exact result. This finding is explained 
in the appendix |Al The underestimation is strongest for 
intermediate interactions and vanishes for weak and for 
strong interactions. However, most importantly, our the- 
ory describes the nature of the ground state correctly 
as a quantum mechanical singlet with proper left-right 
correlations, which manifest themselves as strong anti- 
ferromagnetic correlations in the present model. This 
feature is evident from Fig. [3] where we present the spin- 
spin correlation function Sr^w = {^n\SRSji'\^R) for 
the 4-site chain with U/t=l (top) and U/t=8 (bottom) 
at half filling compared to the exact result. Note that 
one observes an increase of the nearest-neighbor correla- 
tion with increasing [/, which is in accordance with the 
general tendency of the Hubbard model towards anti- 
ferromagnetism at half filling. The correlation length is 
well reproduced within our theory, albeit somewhat un- 
derestimated as compared to the exact result. The lat- 
ter point is to be expected, because the one-dimensional 
Hubbard model is a Luttinger liquid with a slow alge- 
braic decay of spin-spin correlations, while for our local 
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Spin-spin correlation function Sr^r/ 



(^rISrSji'I^r) for the 4-site chain with U/t=l (top) 
and U/t=8 (bottom) at half filling. The many-particle 
wave functions |$h) is that of the Anderson model with the 
impurity at the site of the first site index of Sr,ri . Site one is 
the left terminal site and site two is the left central site. The 
results from our theory, shown as black bars are compared 
with the exact result for the Hubbard chain, shown as white 
bars. The results show the antiferromagnetic coupling of the 
approximate ground states. The antiferromagnetic correla- 
tion grows with increasing U. The correlation length of the 
is well described by our theory, even though underestimated 
for large U. 



approximation we can expect that the corresponding cor- 
relation function is rather that of a Fermi liquid with the 
standard Friedel oscillation exponent. 

A correct description of the ground state, in particular 
reproducing the proper anti-ferromagnetic correlations, 
is by no means trivial. The reason is that left-right 
correlation in a singlet cannot be captured by a single 
Slater determinant. Thus, it cannot be obtained in an 
independent-electron picture. 

In order to compare our results with the independent- 
electron picture, we performed calculations based on a 
static mean-field theory. To generate a static mean-field 



description within our density-matrix functional frame- 
work, one simply has to replace F^'^ [p] in Eq. |43] by the 
expression W^'^'ip] = UpR-f^R-fPm^Ri. 

The ground state in the mean-field approximation 
is a single Slater determinant. For small interaction 
strengths, the wave function is formed from the two 
bonding orbitals and the electrons remain uncorrelated. 
As a consequence the total energy raises with constant 
slope with the interaction strength. At approximately 
U = 2t, the mean-field ground-state wave function under- 
goes a transition to an anti-ferromagnetic state. As in a 
typical second order transition, the magnetization grows 
with an approximate square-root behavior away from the 
transition. In the largc-U limit, each site has the mag- 
netic moment of one electron. The anti-ferromagnetic 
state is a so-called broken-symmetry state, which is not 
an eigenstate of the total spin. For the dimer, such a 
broken-symmetry state can be described as a superpo- 
sition of a singlet and a triplet state. The true ground 
state, however, is a pure singlet, because only the singlet 
state can profit from some remaining covalency. 

Another interesting quantity for the Hubbard model is 
the double occupancy, defined as da = ($|nfl^-|-7ifl^4.|$). 
In the Hubbard model, the interaction energy depends 
exclusively on the sum of the double occupancies. 

In our theory, we extract the double occupancy of a 
given site from the Anderson impurity model with the 
interaction at that site. Hence, the double occupancy is 
calculated for each site from a different Anderson impu- 
rity model and thus from a different wave function. 
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FIG. 4: Double occupancy the 6-site Hubbard chain at half- 
filling as function of the interaction strength from our theory 
(symbols) in comparison with the exact result (lines). Finite- 
size effects can be estimated by comparing the central sites 
(black squares) with the subsurface (open squares) and sur- 
face sites (open circles). The exact result does not distinguish 
appreciably between central and subsurface sites, which indi- 
cates that the underestimation by our theory is not a finite 
size effect. 




FIG. 5: Natural orbitals of the 6-site Hubbard chain for U=0 
(left) and U/t=8 (right). The radius of each sphere is propor- 
tional to the amplitude of the natural orbital, and the area 
is proportional to the site occupation. Full and open circles 
represent the different signs of the orbital coefficients. The 
orbitals are arranged with decreasing occupation from bot- 
tom to top. For each spatial orbital shown, there are two spin 
orbitals. 



The double occupancy for the 6-site Hubbard chain 
is shown in figure 31 The continuous suppression of the 
double occupancy with increasing interaction strength is 
well reproduced. Also the difference between surface sites 
and the center of the chain is in agreement with the exact 
calculation. However, our theory tends to underestimate 
the double occupancy resulting in a somewhat steeper 
slope in the small-U limit. The large-U limit is well re- 
produced. 

The double occupancy provides clear evidence for the 
strength of the presented theory as compared to the 
mean-field theory. In the mean-field theory, the dou- 
ble occupancy remains at d = ^ until U « 2i, from 
where it starts to decrease as the system builds up anti- 
ferromagnetic correlations. In the exact result, and in our 
approximation, the double occupancy starts to decrease 
as soon as the interaction is switched on. The failure of 
mean-field theory is due to its inability to produce the 
correct singlet ground state. 

While the Anderson impurity models provide a local 
perspective and insight on electron correlations as they 
are present in the many-particle wave function, the den- 
sity matrix provides the global perspective. The density 
matrix is fully described by a set of one-particle orbitals, 
the natural orbitals, and their occupations. While the 
natural orbitals are distinct from Kohn-Sham orbitals of 
density-functional theory, they provide similar type of 
information on the system. 

It is striking that the natural orbitals depend very little 
on the interaction strength as seen in figure El Qualita- 
tively there is no difference and the quantitative changes 
are barely visible in this representation. However, the in- 
teraction strength has a strong effect on the occupations 
as shown in figure [SI While the occupations in the (non- 
degenerate) ground state of a non-interacting system are 
always either zero or one, all occupations become frac- 
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FIG. 6: Occupations /„ of the natural orbitals as function of 
the interaction strength from our theory (circles) compared 
to exact result. The occupations become more equally dis- 
tributed with increasing interaction strength. 



tional if the interaction is switched on. From a chemical 
point of view, the depopulation of bonding orbitals and 
the population of formerly unoccupied antibonding or- 
bitals reflects the depression of the covalent interaction 
by the interaction. From a solid-state point of view, the 
behavior seen in Fig. [6] is reminiscent of the momentum 
distribution functions of a Fermi-liquid or the Luttinger 
liquid. While for noninteracting particles the momentum 
distribution has a sharp drop of the occupations from one 
to zero at the Fermi momentum, for the interacting elec- 
tron systemm, it exhibits a gradual decrease with either 
a drop of size Z~^ < 1, the quasi-particle renormaliza- 
tion factor of a Fermi liquid, or, for a Luttinger liquid, a 
sharp change with an algebraic singularity at the Fermi 
momentum 

The smearing out of the occupations implies that the 
number of one-particle wave functions to be considered 
must be larger than that required for density-functional 
theory. If we venture into an extrapolation of our find- 
ings to real systems, we expect that it will be necessary to 
include one-particle wave functions |'0n) up to the anti- 
bonding orbitals of bonds for which left-right correlation 
is important. This implies that we need to include or- 
bitals up to several electron volts above the Fermi-level. 
Working in our favor is that in the strongest bonds, with 
high-lying antibonding states, the kinetic energy is large 
and therefore dominating over the Coulomb repulsion. 
Correlations are dominant in relatively weak bonds, ei- 
ther when they are stretched to the dissociation limit, or 
when they involve d- and f-electron orbitals or n bonds. 

One of the major problems plaguing density-functional 
calculations is the band-gap problem. Here, we do not 
refer to the optical band gap, which is an excited state 
property and which, in a strict sense, is neither accessible 
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FIG. 7: Chemical potentials fi versus filling N/L from our the- 
ory (spheres and dashed lines) and the exact Hubbard chain 
(solid lines) for the four-site (left) and the six-site (right) Hub- 
bard chains as function of the particle number for U/t — 1/2, 
U/t = 3 and U/t — 90 from top to bottom. The dashed lines 
are a guide to the eye and connect the data points. (Color 
online) 



by ground-state DFT nor by our theory. Instead we in- 
vestigate here the thermodynamic gap, which we define 
as 



A = hm 

s^Q \ dN 



N+S 



dE 
dN 



N-S 



Because E{N) is a series of line-segments interpolating 
between the values for integer particle numbers N, this 
definition is identical to the more common expression 
A = £:(iV+l)-2^(iV)+£'(7V-l), which is the difference 
A ~ A — I between electron affinity A and ionization po- 
tential /. Because E{N) is a true ground-state property, 
it should be properly described both by our theory and 
by correct density-functional theory. A failure of describ- 
ing the thermodynamic gap has a fundamental impact on 
the defect chemistry in semiconductors and insulators. 

The chemical potential n{N) = ^ for different in- 
teraction strengths calculated with our theory is shown 
in figure [7] and compared to the exact thermodynamic 
band gap, included as horizontal lines in the figure. It 
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FIG. 8: Total energy for fractional particle numbers N divided 
by the number of sites L of the Hubbard chains for U /t — 1/2 
(first panel), U/t = 3 (second) and U/t = 90 (third). Symbols 
denote results from our theory, lines the exact behavior. 



is obvious that our theory does not yield a discontinu- 
ity. Despite of this, the value of the band band gap can 
nevertheless be obtained using Slater's transition rule as 
A = fj,{N + 5) — fJ-{N — i) from the difference of the 
chemical potentials at half-integer occupations. 

The failure to describe the discontinuity of the chem- 
ical potential, and thus the band gap, reminds of the 
behavior of approximate density-matrix functionals;^ 
These functionals do not make the local approximation, 
but use approximate numerical representations of the ex- 
change correlation energy. Also here, the discontinu- 
ity vanishes, while the band gap calculated with half- 
occupied orbitals are satisfactory predictions. 

Except for the non-interacting case, the total energy 
appears to be a continuous function of the particle num- 
ber within our theory, while the exact result has jumps in 
the slope at integer fillings. We attribute this deficiency 
of our theory to its tendency to underestimate the total 
energy. As seen in figure [8] this underestimation is most 
prominent for fractional particle numbers at intermediate 
interaction strengths. For small or huge values of U the 
effect is less visible; but here, too, the variation across in- 
teger fillings is smooth, as it is apparent from Fig. [71 As 
discussed in appendix [X] this defect is related to the fact 



that the wave function for each Anderson-impurity model 
is optimized individually. An approximation based on a 
common wave function for all sites may direct towards a 
remedy of this problem. 



V. CONCLUSION 

We proposed a pathway to include explicit correlations 
into density-functional based calculations. The guiding 
ideas have been to develop a strict variational approx- 
imation of the interacting electron gas, that avoids the 
dimensional bottleneck of may-particle quantum theory, 
and that allows seamless incorporation into electronic 
structure methodology of density-functional theory. 

Let us summarize what we have accomplished so far: 
Based on a general formulation using a density-matrix- 
functional description of the solid, the many-particle 
problem for the lattice with the full interaction has been 
replaced by a collection of generalized multi-orbital An- 
derson impurity models, defined through the functionals 
F'^'', see Eq. ([26]) . These Anderson impurity models are 
defined relative to the DFT functional as accurate way 
to obtain the quasi-particle dispersion including effects 
of non-local correlations and proper screening. Since 
the local DFT contributions can be explicitly subtracted 
from the density-matrix functional, the problem of dou- 
ble counting can be treated without further assumptions 
in this approach. Restricting the actual correlations to 
one site only has proven to be a powerful idea in many- 
particle schemes such as Dynamical Mean-Field theory. 
We believe that the density-matrix functional approach 
plus local approximation as proposed in this paper, when 
properly implemented, will provide a tool to treat corre- 
lation effects in a spirit similar to DMFT, however using 
the information of the density-matrix only. 

As test case we applied our method to finite Hubbard 
chains as simple model systems. The extremalization of 
the the density-matrix functional is done on the fly using 
an explicit constrained-search algorithnt^. This itera- 
tive approach has been solved by a dynamical approach 
inspired by the Car-Parrinello method. The iterative na- 
ture of the constrained-search algorithm holds promise 
to important cost reductions in an outer self-consistency 
loop. 

In contrast to static mean-field theory, our method pro- 
duces at half filling the correct singlet state with strong 
anti- ferromagnetic correlations, and thus gives a quali- 
tative description of left-right correlation. The local ap- 
proximation smears out the band-gap, which has been 
observed similarly in approximate density-matrix func- 
tionals. The energy levels of the many-particle system 
are, however, well approximated using the chemical po- 
tentials half-integer particle numbers. 

Presently, each of these Anderson impurity models still 
represents a nontrivial many-particle problem on the infi- 
nite lattice. Thus, apparently nothing has been gained so 
far regarding a calculation of physical properties in this 
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approximation. The important next step now is to refor- 
mulate the effective quantum impurity problem in a way 
that it can be solved efficiently. In standard DMFT the 
important aspect is that the actual lattice structure only 
enters through the density of states. As consequence, the 
solution of the effective quantum impurity problem can 
be done without ever resorting to the actual lattice. 

It is obviously important to develop a similar ap- 
proach for the local density-matrix functional. The most 
straightforward idea here would be to treat the quantum 
impurity as open quantum system and apply techniques 
developed in this field to facilitate the creation of an ef- 
ficient impurity solver for our theory. 
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Appendix A: Local approximation 



Let I'I'oip]) be the ground-state wave function for the 
interaction ^^ Wr and a one-particle term Hq so that 



ho[p] + Y,WR-£o[p] 



\Mp]) = 0- 



(A2) 



The one- particle term ft.o[p] is determined such that 
|$o[p]) produces the specified reduced density matrix p. 
Let furthermore |$fl[p]) be the ground-state wave func- 
tion with only a local interaction on site R and with the 
specified reduced density matrix: 



hR[p] + Wr~ £r[p] \<^r[p]) = 



(A3) 



With these definitions, we obtain the desired inequality 
Eq.EIl 



F^'^'^-[p]=J2{Mp] 



Wr 



Mp] 



> Y^UrIp] Wr $fl[p] ) =J2f'^-[p] . (A4) 



Here we show that distributing the interaction into 
many Anderson-impurity models leads to an underesti- 
mation of the energy, i.e. 



F^H Wk [p] > Y^ pWn. [p] 



(Al) 



The argument rests on the simple fact that the energy 
decreases, if one optimizes wave functions for each im- 
purity model individually, rather than to select one wave 
function that needs to suit all local interactions simulta- 
neously. 
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